Resistive transition in frustrated Josephson-junction arrays on a honeycomb lattice 
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We use driven Monte Carlo dynamics to study the resistive behavior of superconducting Joseph- 
son junction arrays on a honeycomb lattice in a magnetic field corresponding to / flux quantum 
per plaquette. While for / = 1/3 the onset of zero resistance is found at nonzero temperature, 
for / = 1/2 the results are consistent with a transition scenario where the critical temperature 
vanishes and the linear resistivity shows thermally activated behavior. We determine the thermal 
critical exponent of the zero-temperature transition for / = 1/2, from a dynamic scaling analy- 
sis of the nonlinear resistivity. The resistive behavior agrees with recent results obtained for the 
phase-coherence transition from correlation length calculations and with experimental observations 
on ultra-thin superconducting films with a triangular pattern of nanoholes. 
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I. INTRODUCTION 

Josephson-junction (JJ) arrays have remarkable prop- 
erties in a magnetic field, which are strongly dependent 
on the geometry of the structure. In addition to be- 
ing realized as two-dimensional arrays of weakly coupled 
superconducting grains^—, they provide important mod- 
els for superconducting wire networks^— and other in- 
homogeneous superconduting systems, when phase fluc- 
tuations of the superconducting order parameter play a 
major role^. An idealized J J array is equivalent to the 
frustrated XY model^, where frustration can be tuned 
by the applied external magnetic field. The frustration 
parameter /, corresponding to the number of flux quan- 
tum per plaquette of the array, sets the average density 
of vortices in the lattice of pinning sites formed by the 
plaquette centers. Depending on the topology of the lat- 
tice of pinning sites and the value of /, a commensu- 
rate vortex lattice is favored in the ground state, allow- 
ing for a phase-coherence transition at finite tempera- 
ture. In this case, the equilibrium phase transitions and 
resistive behavior of the superconducting array are rea- 
sonably well understood for simple low-order commen- 
surate phases such as / = 1/2 on a square array^ and 
/ = 1/3 on a honeycomb arrayi -. The magnetoresis- 
tance for a square JJ array, for example, oscillates with 
the applied magnetic fieldi^, displaying minima at in- 
teger values of / and secondary minima at / = 1/2 for 
decreasing temperatures, corresponding to resistive tran- 
sitions at different temperatures^. The onset of zero- 
resistance for decreasing temperatures marks the phase- 
coherence transition in the JJ array, which for integer 
/ is expected to be in the Koterlitz-Thouless (KT) uni- 
versality class. Dynamical transitions under an external 
driving current have also been studied for / = 1/2 on 
a square lattice, leading to interesting nonequilibrium 
phase diagrams 11 . However, when the vortex lattice is 
incommensurate with the pinning sites, as for irrational 
/ on a square J J array2^&i2~— or / = 1/2 on a hon- 
eycomb JJ arra y 7 ' 10 ' 17 " —, the possible phase transitions 
are much less understood, showing some features of a vor- 



tex glass without disorder and dynamical freezing at low 
temperatures. In particular, a JJ array on a honeycomb 
lattice with / = 1/2, should display interesting resistive 
behavior. As a model of phase fluctuations, it should be 
relevant to ultra-thin superconducting films with a peri- 
odic pattern of nanohole a 20 ' 21 , which can be regarded as 
a lattice of pinning centers. While for a square lattice 
of nanoholes, the magnetoresistance oscillates with the 
applied field, displaying secondary minima at / = 1/2 as 
for a square JJ array™ for a triangular lattice 2 ^ it shows 
only minima at integer flux quantum per lattice unit cell. 

In early Monte Carlo (MC) simulations of the fully 
frustrated XY model on a honeycomb lattice^, a phase- 
coherence transition at a nonzero temperature in the KT 
universality class was suggested and therefore a resistive 
transition would be expected for a JJ array in the same 
lattice with / = 1/2. On the other hand, a different 
calculation 17 suggested a spin-glass like transition. It 
was also suggested^ that only a crossover region rather 
than an equilibrium phase transition should occur at any 
nonzero temperature. Recently^, it was argued that 
vortex-ordered phases could be possible at nonzero tem- 
peratures but for very large systems, beyond the ones cur- 
rently studied numerically or even experimentally. How- 
ever, the question of the resistive transition was not inves- 
tigated. In a recent MC study of phase coherence in the 
fully frustrated XY model a zero-temperature transition 
scenario^ was proposed, where T c = but the divergent 
correlation length, £ oc T _I/ , should lead to measurable 
effects at finite temperatures in the linear and nonlinear 
resistivity, determined by the thermal critical exponent 
v. So far, a direct calculation of the resistive behavior 
of JJ arrays on a honeycomb lattice and comparison to 
experiments have not been presented. 

In this work, we present results for the resistive behav- 
ior obtained by driven Monte Carlo dynamics. While for 
/ = 1/3 a resistive transition is found at nonzero temper- 
ature, for / = 1/2 the results are consistent with a transi- 
tion scenario where the critical temperature vanishes and 
the linear resistivity shows thermally activated behavior. 
We determine the thermal critical exponent v of the zero- 
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FIG. 1: JJ array on a honeycomb lattice. Filled circles 
represent superconducting grains and the lines the Josephson 
junctions between them. 

temperature transition for / = 1/2, from a dynamic scal- 
ing analysis of the nonlinear resistivity. Its value is in fair 
agreement with recent calculations for the frustrated XY 
model from finite-size correlation length scaling^. A dy- 
namical freezing at lower temperatures is also identified 
from deviations of the fluctuation-dissipation relation be- 
tween linear resistivity and voltage autocorrelations. The 
resistive behavior is consistent with some experimental 
observations in ultra-thin superconducting films with a 
triangular lattice of nanoholes^i, taking into account the 
effects of weak Josephson-coupling disorder. 

II. MODEL AND DRIVEN MONTE CARLO 
SIMULATION 

We consider a JJ array in a uniform transverse mag- 
netic field described by the Hamiltonian 

H = - J H c°s(0* " Oj - Mi)- Jj2i0i - i+i ), (1) 

<ij> i 

where 6*; is the phase of the local superconducting or- 
der parameter of the grains located on the sites of a 
two-dimensional honeycomb lattice with lattice spacing 
a, as illustrated in Fig. [1] The first term is the con- 
tribution from the Josephson-coupling energy between 
nearest neighbor grains. For uniform coupling we set 
Jij = Jo, a constant independent of the magnetic field. 
The line integral of the vector potential Ay due to the 

external field B — V x A is constrained to Aij = 
around each hexagonal plaquette, where / is the number 
of flux quantum <p = hc/2e per plaquette. This model 
is periodic in / with period / = 1. In the calculations 
we choose a gauge where Ay = 2irfrii/2 on the (tilted) 
bonds along the horizontal rows numbered by the inte- 
ger m and Ay = on the vertical bonds of the lattice. 
The second term in Eq. ([lj represents the effects of an 



o.i 



0.01 - 
0.001 - 



-»-T= 


0.26 




0.25 


-*-T= 


0.24 




0.235 




0.23 


-e-T= 


0.2275 


-B-T= 


0.226 


^T= 


0.225 


^rT= 


0.224 


^T= 


0.223 


-»-T= 


0.2225 



io- 5 t LI , , , J 

0.001 0.005 0.010 0.050 0.100 0.500 

J 

FIG. 2: Nonlinear resistivity E/J as a function of current 
density J at different temperatures T, for / = 1/3. System 
size L = 60. The dashed line indicates power-law behavior 
E oc J 3 . 

external driving current density (2e/h)J applied in the x 
(horizontal) direction, coupling to the phase difference, 
9i — 8i + x, between nearest neighbors sites in this direc- 
tion. When J / 0, the total energy is unbounded and the 
system is out of equilibrium. The lower-energy minima 
occur at phase differences Qi — 9i +x which increases with 
time t, leading to a net phase slippage rate proportional 
to < d{9i — 9i + x)/dt >, corresponding to the voltage 
Vij+s. For convenience, we use units where 2e/h — 1, 
J a = 1 and o = l. 

To study the current- voltage behavior, we use a driven 
MC dynamics method^. The time dependence is ob- 
tained by identifying the MC time as the real time t and 
we set the unit of time dt = 1, corresponding to a com- 
plete MC pass through the lattice. For convenience, the 
honeycomb lattice is defined on a rectangular geometry 
(Fig. [T]), with linear size given by a dimensionless length 
L. In terms of L, the linear size in the x and y directions 
can be written as L x — L\fia and L y = to, respectively. 
This corresponds to 2L junctions along the horizontal 
rows. The usual periodic boundary conditions are used 
in the y-direction and periodic (fluctuating twist) bound- 
ary conditions^ in the x direction. The twist boundary 
condition adds a new dynamical variables u x , correspond- 
ing to a uniform phase twist between nearest-neighbor 
sites along the x direction. A MC step consists of an at- 
tempt to change the local phases 9i and the phase twist 
u x , using the Metropolis algorithm. If the change in en- 
ergy is AiT, the trial move is accepted with probability 
min{l, exp(— AH/kT)}. The external current density J 
in Eq. ^ biases these changes, leading to a net voltage 
(phase slippage rate) across the system in the s-direction 
given by 

V = 2Lj t u x , (2) 
in arbitrary units. Compared to the usual Langevin 
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III. RESULTS AND DISCUSSION 
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FIG. 3: Nonlinear resistivity E/J as a function of current 
density J at different temperatures T, for / = 1/2. System 
size L — 60. 



dynamics^, this MC method allows to access much 
longer time scales, which is required to obtain reliable 
data at lower temperatures and current densities. We 
have determined the electric field E = V/(2L) and non- 
linear resistivity p = E / J as a function of the driving 
current density J, in the x direction, for different tem- 
peratures T and different system sizes L. We used typi- 
cally 5 x 10 6 MC steps to reach the nonequilibrium steady 
state and equal time steps to perform time averages, with 
additional averages over 6—12 independent runs. In a 
MC step, the maximum changes in the local phases 0, 
and the phase twist u x were fixed to ±7r and ±7r/(2£), 
respectively. 

The linear resistivity, pl = lim,/_>o E/ J, can be deter- 
mined from the nonlinear behavior p(J), obtained from 
the driven MC simulations, by extrapolating the numeri- 
cal results to vanishing currents. It can also be obtained, 
independently, from equilibrium voltage fluctuations and 
therefore can be calculated in absence of an imposing 
driving current (J = 0). From Kubo formula, the lin- 
ear resistance is given in terms of the equilibrium voltage 
autocorrelation as 



(3) 



Since the total voltage V is related to the phase difference 
across the system A6(t) by V — dA9(t)/dt, we find more 
convenient to determine Rl from the long-time equilib- 
rium fluctuations^ of A8(t) as 



R L = —((A9(t)-A0(O)f), 



which is valid for sufficiently long times t. 



(4) 



First, we consider the resistive behavior when / = 1/3. 
For this value of the frustration, it is known that a hexag- 
onal vortex lattice commensurate with the honeycomb 
lattice is the ground stated, and therefore a resistive 
transition would be expected at a temperature smaller 
or equal the vortex lattice melting. Fig. [2] shows the 
nonlinear resistivity E/J as a function of temperature, 
for the largest system size L = 60, where finite-size ef- 
fects are small. For decreasing current densities J, the 
nonlinear resistivity E/J tends to a finite value at high 
temperatures, corresponding to the linear resistivity pL, 
but extrapolates to very low values at lower tempera- 
tures. This behavior is consistent with a resistive tran- 
sition occurring at a critical temperature in the range 
T c (f = 1/3) = 0.224-0.225. In fact, it is slightly smaller 
than the vortex lattice melting transition estimated from 
recent equilibrium MC simulations of the frustrated XY 
model on a honeycomb lattice^, T m — 0.226(1). At the 
resistive transition, a power-law relation E cx J z+1 is 
expected at sufficiently small currents from the scaling 
theory----, where z is the dynamical critical exponent. For 
the usual KT transition it is known2i2& that 2 = 2. In 
the present shown by the dashed line in Fig. [3J a 

power-law separating the T > T c from T < T c behavior 
at small currents is compatible with = 2. However, 
further work taking into account finite-size effects is re- 
quired to investigate the critical behavior in detail. In 
any case, the above results show clear evidence of a re- 
sistive transition at finite temperature for / = 1/3. 

In contrast to the resistive behavior in Fig. [21 when 
/ = 1/2 the nonlinear resistivity E/J tends to a finite 
value for decreasing currents even at low temperatures as 
shown in Fig. [3] Although we can not exclude a tran- 
sition at much lower temperatures, where reliable data 
could not be obtained as discussed below, this behavior 
is consistent with a resistive transition occurring only at 
zero temperature. Recent equilibrium MC simulations 
suggested such zero-temperature transition scenario^, 
where T c = for the phase-coherence transition but the 
finite correlation length for T > leads to measurable 
effects in the nonlinear resistivity. In fact, the behav- 
ior in Fig. [3] has the main features expected for a zero- 
temperature resistive transition. The linear resistivity 
Pl, corresponding to zero current limit of E/J, decreases 
rapidly with decreasing temperature and for increasing J, 
E/J cross over to a nonlinear behavior at a characteristic 
current density J„/ , which also decreases with decreasing 
temperature. 

To verify in which temperature range the values ap- 
proached at low currents in Fig. [3] correspond indeed to 
the linear resistivity pl, we show in Fig. 3] the tem- 
perature dependence of pl obtained from the nonlin- 
ear resistivity as pl — lirnj_>o-E/J and, without cur- 
rent bias, from Eq. ([3]). These values obtained from 
nonequilibrium and equilibrium calculations agree with 
each other above a temperature Tf ~ 0.11 and deviate 
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FIG. 4: Temperature dependence of the linear resistivity pL 
for / = 1/2 obtained from nonlinear resistivity (J — > 0) and 
from voltage fluctuations (J = 0). System size L = 60. The 
separation of the curves gives an estimate of the dynamical 
freezing temperature Tf. The dashed line is an Arrhenius fit 
for T > Tf. 
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FIG. 6: Scaling plot of the nonlinear resistivity E/J from 
Fig. [3l system size L = 60, for different temperatures 
above Tf and current range J < 0.03, giving the estimate 
v = 1.40(9). Inset: Finite-size scaling plot of the crossover 
current J n i for different temperatures T and system sizes L, 
giving the estimate v = 1.15(9). The dashed line is the fit 
used in the data collapse procedure. 
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FIG. 5: Temperature dependence of the crossover current 
J„i for / = 1/2 and / = 1/3. System size L = 60. The 
dashed line is a power-law fit to J n ; oc T 1+ " , giving the esti- 
mate v = 1.17(14). The arrow indicates the estimated critical 
temperature for / = 1/3. 

significantly at lower temperatures. Since this agreement 
is only expected when the voltage autocorrelation in Eq. 
(j3]) is obtained in true equilibrium, one can regard Tf 
as a signature of a dynamical freezing transition, below 
which equilibrium is not achieved due to very large re- 
laxation time. Interestingly, a dynamical freezing tran- 
sition near the same temperature was also identified in 
recent equilibrium MC simulations by other methods and 
different dynamics^. The apparent KT transition^ and 
spin-glass transition 1 - observed in earlier MC simulations 
could be attributed to slow dynamics effects of such dy- 
namical freezing. 

The straight-line behavior of Pl{T) for T > Tf in 



the log-linear plot of Fig. 2] indicates an activated Ar- 
rhenius behavior, where the linear resistivity decreases 
exponentially with the inverse of temperature with a 
temperature-independent energy barrier, estimated as 
Eb = 1.16(4)J Q . If such behavior extrapolates to lower 
temperatures, it suggests that the linear resistivity can be 
very small but nevertheless remains finite for decreasing 
temperatures and therefore there is no resistive transition 
at finite temperatures. However, as will be described 
below, the system behaves as if a resistive transition 
occurs at zero temperature, corresponding to a phase- 
coherence transition where the critical temperature van- 
ishes, T c = 0. 

A detailed scaling theory^ of the resistive transi- 
tion with T c = has been described in the con- 
text of the current- voltage characteristics of vortex-glass 
models^— of disordered two-dimensional superconduc- 
tors but the arguments should also apply to the present 
case. The basic assumption is the existence of a second- 
order phase transition. The correlation length £ is finite 
for T > but it increases with decreasing temperature 
as £ oc T~ v , with v a critical exponent. The divergent 
correlation length and relaxation time r near the transi- 
tion determine both the linear an nonlinear resistivity 
behavior leading to current-voltage scaling sufficiently 
close to the critical temperature and sufficiently small 
driving current. If the data satisfy such scaling behav- 
ior for different driving currents and temperatures, the 
critical temperature and critical exponents of the under- 
lying equilibrium transition at J = can then be deter- 
mined from the best data collapse. The dimensionless 
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FIG. 7: Scaling plot of the nonlinear resistivity El J as in 
Fig. [6] for a larger system size L = 72, giving the estimate 
v = 1.36(8). Inset: Same scaling plot for L = 90, giving the 
estimate v = 1.33(8). 



ratio E/JpL should satisfy the scaling forrr>2& 
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(5) 



where g is a scaling function with g(0) = I. A crossover 
from linear behavior, when g(x) ~ 1, to nonlinear be- 
havior, when <?(:c) >> 1, occurs when x ~ 1 which leads 
to a crossover current density at which nonlinear behav- 
ior sets in, decreasing with temperatures as a power law, 
Jul oc T/i oc T 1+v . The scaling form in Eq. © con- 
tains a single critical exponent v and does not depend 
on the particular form assumed for the divergence of the 
relaxation time r. However, for sufficiently low tempera- 
tures, the relaxation process is expected to be thermally 
activated^ with r oc exp(EbfkT). This corresponds for- 
mally to a dynamic exponent z — > oo, if power-law be- 
havior is assumed for the relaxation time t a f . The 
linear resistivity should scale as 2 ^ pl oc 1/r and there- 
fore it is also expected to have an activated behavior, 
PL oc exp(— Eb/kT). In general, the energy barrier Eb 
also scales with the correlation length as Eb oc £^ , which 
leads to a temperature-dependent barrier Eb oc T~^ v . A 
pure Arrhenius behavior corresponds to if) = 0. 

The behavior of the nonlinear and linear resistivity in 
Figs [3] and [4] above the dynamical freezing temperature 
Tf are quite consistent with the predictions from the scal- 
ing theory. Fig. [5] shows the temperature dependence 
of the crossover current J n i, defined as the value of J 
where Ej Jp^ starts to deviate from a fixed value, cho- 
sen to be c = 1.2. For the lowest temperature range 
above Tf, the linear behavior in the log- log plot is con- 
sistent with the expected power-law J„; oc T 1+I/ for a 
zero-temperature transition. From the power-law fit we 
obtain a first estimate of the exponent v = 1.17(14). In 
contrast, for / = 1/3, the behavior in the lowest tem- 
perature range does not allow a similar power-law fit; 



J n i curves down for decreasing temperatures and extrap- 
olates to zero at a finite temperature, consistent with 
a resistive transition at a nonzero critical temperature 
found for / = 1/3. The nonlinear resistivity data also 
satisfies the scaling form for different driving currents and 
temperatures. Fig. [5] shows a scaling plot of the nonlin- 
ear resistivity above Tf according to Eq. ((5]), for a large 
system size L = 60, where the finite-size dependence is 
small. The best data collapse provides an estimate of 
the critical exponent v = 1.40(9). The data collapse is 
achieved quantitatively by means of a least-squares fit 
metho d 19 ' 28 , varying the parameter v. The scaling func- 
tion g{x) is approximated by a Taylor series expansion for 
small x, truncated beyond 4th order, which is used to fit 
the data and provide the least-square residuals. The er- 
ror estimate here corresponds to the statistical error from 
the least-squares method and do not include systematic 
effects. To check for systematic errors from finite-size ef- 
fects, which were assumed negligible in the scaling form of 
Eq. ([5]), the same data collapse procedure was repeated 
for larger system sizes, as shown in Fig. [7J The results 
of these estimates, v = 1.36(8) for L = 72 and 1.33(8) 
for L = 90, agree within the statistical errors but indi- 
cate that the central estimate of v decreases slowly with 
system size. The nonlinear resistivity should also satisfy 
the expected finite-size behavior in smaller system sizes 
when the correlation length £ approaches the system size 
L. According to finite-size scaling, the scaling function 
in Eq. ^ , should also depend on the dimensionless ratio 
L/£ and so, to account for finite-size effects, the nonlinear 
resistivity should satisfy the scaling form 



E 
J PL 



J 



,L 1 ' U T). 



(6) 



The scaling analysis of the whole nonlinear resistivity 
data is rather complicated in this case since the scal- 
ing function depends on two variables. To simplify the 
analysis 2 ^ we first estimate the temperature and finite- 
size behavior of the crossover current density J n i where 
nonlinear behavior sets in, as the value of J where 
E/JpL = c, a constant. Then, from Eq. ©, the finite- 
size behavior of J n \ can be expressed in the scaling form 



J nl L^' v = =g{L 1 ' v T). 



(7) 



The best data collapse according to the scaling in Eq. 
([7]) provides an independent estimate of the critical ex- 
ponent v. The inset in Fig. [5] shows that indeed the 
values of J n i for different system sizes and temperatures 
satisfy this scaling form with v = 1.15(9). To check for 
systematic errors due to corrections to finite-size scaling, 
the data collapse was repeated dropping the smaller sys- 
tem sizes. Dropping system size L — 24 gives v = 1.17(7) 
and L = 24 to L = 36 gives v = 1.16(7). Since the re- 
sulting changes are small compared with the errorbars, 
systematic errors of this kind are not significant for this 
range of system sizes. The two independent estimates 
of v obtained above, 1.33(8) from the largest system 
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FIG. 8: Temperature dependence of the linear resistivity p^ 
for different frustration parameters / and coupling disorder 
strengths D. System size L = 60. 



size and 1.15(9) from finite-size scaling, are not compat- 
ible within the estimated errors. However, the former 
value could still be affected by finite-size effects. The 
latter value should be more accurate since it is based 
on finite-size scaling. This value is in reasonable agree- 
ment, within the estimated errors, with the critical ex- 
ponent for the zero-temperature phase-coherence tran- 
sition, v p h — 1.29(15), of the frustrated XY model ob- 
tained recently by correlation length calculations using 
equilibrium MC simulations^. The Arrhenius behavior 
for the linear resistivity pl in Fig. 0] is also consistent 
with the exponential divergence of the relaxation time r 
found in the equilibrium MC simulations. 

Some experimental observations on ultra-thin super- 
conducting films with a triangular pattern of nanoholes 21 
are consistent with the zero-temperature resistive transi- 
tion for / = 1/2. In the regime where phase fluctuations 
of the superconducting order parameter are more impor- 
tant than amplitude fluctuations^, this system can be 
described by an array of superconducting "grains" cou- 
pled by Josephson junctions in a suitable geometry. The 
simplest model consists of a Josephson-j unction array 
on a honeycomb lattice, with the triangular lattice of 
nanoholes corresponding to the lattice of pinning sites 
(plaquette centers in Fig. [IJ and the number of flux 
quantum per unit cell of the nanohole lattice correspond- 
ing to the frustration parameter / of the array. In fact, 
the measured resistance of samples which are supercon- 
ducting at low temperatures and low magnetic fields, os- 
cillates as a function of the magnetic field, displaying 
minima at integer values of / but no secondary minima 
at/ = 1/2, as expected from the present results for the 
honeycomb J J array. Moreover, for / = 1/2, the temper- 
ature dependence of the resistance shows the expected 
Arrhenius behavior, consistent with a vanishing critical 
temperature. However, the measured magnetoresistance 
does not display minima at / = 1/3, which would be ex- 
pected from the above calculations for temperatures near 



T c (f — 1/3). Although the available temperatures in the 
experiments may not be sufficiently small to observe this 
feature, it could also be the effect of quenched disorder 
in the Josephson couplings. In fact, it was recently sug- 
gested that inhomogeneities in the film thickness could 
lead to significant variations in the weak links between 
superconducting islands^. 

We have performed additional calculations to verify 
the qualitative effect of weak disorder of the Josephson 
couplings on the magnetoresistive behavior. We consider 
a simple random-coupling model, where Jij in Eq. (JT|) is 
defined as = J (l ± D), with equal probability, and 
disorder strength parameter D. The JJ array is still as- 
sumed to be on a perfect honeycomb lattice. The resistiv- 
ity as a function of temperature was calculated by averag- 
ing over different realizations of the disorder. Fig. [8]com- 
pares the temperature dependence of pl obtained with- 
out current bias, from Eq. ([3]), for / = 0, / = 1/3 and 
/ = 1/2, and different disorder strengths D. While the 
behavior characteristic of a finite-temperature transition 
for / = and zero-temperature transition for / = 1 /2 
remain for increasing disorder, the resistive behavior for 
/ = 1/3 changes to an Arrhenius form above a disorder 
strength D ~ 0.35. In this case, the magnetoresistance 
should only display minima at integer values of /, as 
observed experimentally^, which in turn suggests that 
coupling disorder should also play an important role in 
modeling other phase coherence properties of this sys- 
tem. When comparing the disorder strength D in the 
model with the thickness variations in the experimen- 
tal system, geometrical disorder in the JJ array due to 
spatial irregularities of the system, should also be taken 
into account. Weak positional disorder of the grains, for 
example, has significant effects both on phase-coherence 
and vortex order at nonzero values of frustration^—, 
even when the Josephson coupling is uniform (D = 0). 
One then would expect that the combined effect of ge- 
ometrical and Josephson coupling disorder in the model 
will result in an Arrhenius behavior for f — 1/3 occurring 
at much lower values of D. These interesting effects and 
a more quantitative comparison of such disorder model 
with the experimental system deserves further work. 



IV. CONCLUSIONS 

We have investigated the resistive behavior of Joseph- 
son junction arrays on a honeycomb lattice using driven 
MC dynamics, focusing mainly on the / = 1/2 frustra- 
tion and its relation to experiments on ultra-thin super- 
conducting films2i. For / = 1/3, a resistive transition 
is found at nonzero temperature, as expected from early 
results of equilibrium MC simulations—. The estimated 
critical temperature is slightly below the melting tran- 
sition of the commensurate vortex lattice^, suggesting 
two separated transitions. However, further work is re- 
quired to obtain a more accurate estimate and to investi- 
gate the critical behavior in detail. For / = 1/2, the re- 
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suits are consistent with a transition scenario where the 
critical temperature vanishes and the linear resistivity 
shows thermally activated behavior. The thermal crit- 
ical exponent v of the zero-temperature transition esti- 
mated from a dynamical scaling analysis is in fair agree- 
ment with recent calculations from finite-size correlation 
length scaling^. A dynamical freezing at a lower temper- 
ature Tf was identified from deviations of the fluctuation- 
dissipation relation between linear resistivity and voltage 
autocorrelations. It should be pointed out that, since 
equilibrium data could not be obtained below Tf, a re- 
sistive transition at much lower temperatures can not be 
ruled out. Moreover, since the scaling analysis assumes a 
second-order phase transition, a first-order resistive tran- 
sition near or below Tf is also not excluded. The resis- 
tive behavior is qualitatively consistent with experimen- 
tal observations in ultra-thin superconducting films with 



a triangular lattice of nanoholes? 1 ., taking into account 
the effects of weak Josephson-coupling disorder. A more 
quantitative comparison to the experimental system, in- 
cluding geometrical disorder— ~—, and the relation be- 
tween the resistive behavior and the vortex structure* 8 
for / = 1/2, as well as / = 1/3, require further work. 
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